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The turbulent transport of impurity particles in plasma edge turbulence is investigated. The 
impurities are modeled as a passive fluid advected by the electric and polarization drifts, while the 
ambient plasma turbulence is modeled using the two-dimensional Hasegawa-Wakatani paradigm for 
resistive drift- wave turbulence. The features of the turbulent transport of impurities are investigated 
by numerical simulations using a novel code that applies semi-Lagrangian pseudospectral schemes. 
The diffusive character of the turbulent transport of ideal impurities is demonstrated by relative- 
diffusion analysis of the evolution of impurity puffs. Additional effects appear for inert ial impurities 
as a consequence of compressibility. First, the density of inert ial impurities is found to correlate with 
the vorticity of the electric drift velocity, that is, impurities cluster in vortices of a precise orientation 
determined by the charge of the impurity particles. Second, a radial pinch scaling linearly with the 
mass-charge ratio of the impurities is discovered. Theoretical explanation for these observations is 
obtained by analysis of the model equations. 



I. INTRODUCTION 

One of the persistent problems in the development of a 
magnetic fusion reactor is the degradation of plasma con- 
finement caused by turbulent transport. In spite of the 
strong confining magnetic field, high levels of energy and 
charged-particle transport across magnetic field lines are 
observed in a variety of experimental configurations, such 
as tokamaks and stellar ators. This transport has received 
the adjective "anomalous," as its features cannot be ex- 
plained by the classical or neoclassical collisional theo- 
ries. Instead, anomalous transport is generally attributed 
to small-scale, electrostatic plasma turbulence. ^'^ In par- 
ticular, the transport in the plasma edge is believed to 
be largely controlled by drift- wave turbulence, which ap- 
pears naturally in this region as a result of the strong 
radial pressure gradient. The drift- wave transport mech- 
anism has been the subject of much theoretical and ex- 
perimental research over the past several decades, and 
was recently reviewed by Horton.^ Nevertheless, much of 
this transport phenomenon is yet to be fully understood. 

Another important problem in the quest for controlled 
thermonuclear fusion that is linked to turbulence con- 
cerns the impurity content of the plasma. The properties 
of fusion plasmas are strongly influenced by the presence 
of impurities. For instance, impurities enhance radiative 
energy losses and dilute the hydrogen fuel within the hot 
plasma center, at the same time being the main agents in 
the erosion and deposition processes that affect plasma- 
facing components.^ Consequently, the understanding of 
the transport mechanisms that control the impurity con- 
tent in the plasma is of crucial importance for the per- 
formance and safe operation of future fusion devices. 

The transport of impurity particles across the mag- 
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netic field is related to that of the bulk plasma, since the 
motion perpendicular to the magnetic field of all charged 
particles in the plasma is dominated by one same velocity 
field, the electric drift velocity. In this sense, impurities 
and bulk plasma are partly subject to a common mech- 
anism of turbulent transport, which relates the problem 
of confinement to that of controlling the impurity con- 
tent. However, the distinctive nature of impurity parti- 
cles, notably their mass and electric charge, as well as 
their different influence on the plasma allow for dispari- 
ties between the transport of impurities and that of the 
bulk plasma. In fact, recent experiments not only re- 
veal disparities in the transport of impurities and bulk 
plasma,^ but also in that of different impurity species. ^« 
In such experimental investigations impurity transport is 
usually characterized in terms of effective diffusion coeffi- 
cients and drift velocities, the latter also known as pinch 
velocities. ^'^ These transport coefficients are typically de- 
termined by analysis of the temporal evolution of injected 
impurity puffs. While the explanation for the diffusive 
character of turbulent transport goes back to the work 
of Taylor and Richardson in the 1920s, the origin of 
the anomalous pinch velocities is still a matter of intense 
research.^ Moreover, the inward pinch observed in lab- 
oratory experiments tends to concentrate the impurities 
in the center of the plasma column,^ precisely where they 
are least wanted, and thus constitutes one of the major 
unresolved transport problems in fusion plasma physics. 

In this contribution we investigate the turbulent trans- 
port of impurities in a generic model for plasma edge tur- 
bulence. Specifically, the investigations are based on the 
simple and extensively studied model for resistive drift- 
wave turbulence proposed by Hasegawa and WakataniA 
For the impurities, we use a consistent passive- fluid 
model that takes into account impurity-particle inertia, 
that is to say, it contemplates the polarization drift. 
This enables the model to exhibit different transport fea- 
tures depending on the mass and charge of the impuri- 
ties, as observed in laboratory experiments. In particu- 
lar, the polarization drift brings compressibility into the 
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model, giving rise to subtle yet qualitatively important 
effects. The investigations are carried out by numerical 
simulations of the model equations using a newly devel- 
oped code that applies semi-Lagrangian pseudospectral 
schemes. We proceed separately for ideal (i.e., massless) 
impurities and inert ial ones. In the case of ideal impuri- 
ties, we focus on the phenomenon of turbulent diffusion, 
which is studied from the relative-diffusion perspective. 
For inert ial impurities, we concentrate on the new effects 
that arise from compressibility. First, the density of in- 
ertial impurities is found to correlate with the vorticity 
of the electric drift velocity, that is, impurities cluster in 
vortices of a precise orientation determined by the charge 
of the impurity particles. This behavior is a direct conse- 
quence of compression by the polarization drift and tur- 
bulent mixing. Second, we discover a radial pinch that 
scales linearly with the mass-charge ratio of the impuri- 
ties. An understanding of this phenomenon is obtained 
by analysis of the evolution equation for the global im- 
purity flux. 

This paper is organized as follows. In Sec.^the model 
for resistive drift-wave turbulence and the passive-fluid 
model for the impurities are described. The diffusive ef- 
fect of turbulence is demonstrated in Sec. lIIII for the case 
of ideal impurities. The effects that arise from impurity 
inertia, the clustering and the anomalous radial pinch, 
are respectively presented in Sees. II VI and IVI The results 
are summarized and discussed in Sec. lVIl Lastly, the Ap- 
pendix is devoted to a brief description of the numerical 
schemes applied in the simulations. 

II. MODEL EQUATIONS 

In this section we provide a brief presentation of the 
resistive drift-wave paradigm used for modeling the tur- 
bulence at the plasma edge. We also describe the passive- 
fluid model for the impurity species, which takes into ac- 
count impurity-particle inertia. 

A. Resistive drift-wave turbulence 

As a basic model for plasma edge turbulence we 
consider the two-dimensional (2D) Hasegawa-Wakatani 
(HW) model for resistive drift- wave turbulence:^ 

Dt{n - x) = C{(p - n) + jj.n'^^n, (la) 

BtU = C{(p - n) + /ia;Vicj, (lb) 

where 

CO = Vi(/9, (Ic) 

Dt = St + z X V±(p'V±. (Id) 

Here, n denotes the fluctuating component of the plasma 
density, cp is the electrostatic potential, and co is the vor- 
ticity of the electric drift velocity ve = ^ X V _\_(p. 



The HW model assumes an unperturbed, uniform 
magnetic field in the z-direction, that is, B = Bqz. The 
X- and ^/-directions correspond respectively to the radial 
and the poloidal directions in a slab at the edge of the 
confined plasma. The electrons are assumed isothermal 
with temperature Tg, whereas the ions are considered 
cold. The equilibrium plasma density no is assumed to 
be homogeneous in the poloidal direction and to decrease 
in the radial direction in such a way that the equilib- 
rium density length scale Ln = no/|dno/da::| is constant. 
Hence, in this thin-layer approximation the second term 
on the left-hand side of Eq. (|Ta)) actually represents the 
advection of the inhomogeneous equilibrium plasma den- 
sity by the electric drift. 

The quantities in Eqs. (^) have been made dimension- 
less using the gyro-Bohm normalization: 

— ^x — ^ t—^t — — ^^^^ 

Ps ' ps ' Ln ' no ps ' Te ps 

Here, Cg = {Te/rUiY^'^ is the sound speed and ps = 
{rUiTeY^'^ / cBq is the the drift scale, that is, the ion gyro- 
radius at the sound speed. The parameter C coupling 
Eqs. (|Ta)) and (|Tb|) is known as the adiabaticity parame- 
ter, and is defined as 

~ e^T^.no Cs/Ln' 

Here, 77,1 is the parallel resistivity and /c,, is the single 
parallel wave number that is considered in the reduc- 
tion to a 2D model. In the strongly collisional limit 
C ^ 0, Eq. ([Tb|) reduces to the 2D Navier-Stokes vor- 
ticity equation. In this so-called hydrodynamic limit the 
equations become decoupled and the plasma density is 
passively advected by the electric drift velocity. In the 
weakly collisional limit C ^ 00, the adiabatic electron 
response yields Boltzmann-distributed density perturba- 
tions, n ^ (p, and the HW model reduces to the well- 
known Hasegawa-Mima equation. 

The diffusive terms in Eqs. (|Ta|) and (flbj) correspond 
respectively to the collisional diffusion of electrons and 
the ion viscosity. These dissipative processes are essen- 
tial for the stabilization of resistive drift modes as well 
as for the feasibility of numerically simulating the HW 
model. Linear analysis of the HW model reveals the ex- 
istence of unstable modes for < C < 00, with all modes 
being unstable in the absence of dissipation^AiiS It is the 
combination of this linear instability with nonlinear in- 
teraction through the advective terms and dissipation at 
small scales that leads to quasistationary turbulent states 
within the HW model. In particular, when initialized 
with a broadband set of small-amplitude, random-phase 
modes, the system evolves in two stages: first, a linear 
phase in which the linearly unstable modes grow and the 
damped ones decay; second, a phase of nonlinear satura- 
tion leading to a quasistationary turbulent state whose 
statistics solely depend on the parameters of the HW 
model. 
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Two relevant nonlinear invariants in the HW model 
are the total energy E and the generalized enstrophy [/, 
defined as 



{n-uofdyi. 



The evolution of these quantities within the HW model 
is governed by 



dE 
d[/ 



— J- n ^ 1 



where 



Tc = C J {n- iff dx, 

= J (/in|V^n|' + /i^|V^cj|')dx, 

= J V±{n - co) ' V±{jj.nn - jj^ojoo) dx. 

Hence, is the only source in the system, and corre- 
sponds to energy extracted from the equilibrium density 
gradient. Conversely, Tc is a sink corresponding to re- 
sistive dissipation through the parallel current, while 
and are sinks arising from collisionality and viscosity. 
Being the sole source of energy and generalized enstrophy 
source within the HW model, is in practice positive 
definite. This fact is used in Sec.0to identify a source 
for radial drift of impurities. 



turbulence model, the flow velocity of the impurities is 
hence given by 



where 



me e ps 



qe rui Lr, 



Here, me and qe are respectively the mass and the electric 
charge of the impurity particles. The impurity model 
resulting from the continuity equation is thus given by 



(2) 



where the term on the right-hand side arises from colli- 
sional diffusion. 

The polarization drift, which represents the effect of 
impurity-particle inertia, constitutes a higher-order cor- 
rection to the electric drift velocity. The influence of 
the former on the dynamics of impurities is determined 
by the dimensionless parameter which varies with the 
mass-charge ratio of the impurity particles relative to 
that of the bulk-plasma ions. For ideal impurities C = 
and the impurities are solely advected by the electric 
drift. For inertial impurities has a nonzero value and 
the polarization drift can play a role in the dynamics. 
The value of is in any case small, since the length scale 
quotient ps/Ln is by assumption small. Nevertheless, the 
polarization drift adds important qualitative features to 
the dynamics of impurities. In the case of static, uni- 
form magnetic field the electric drift is incompressible, 
whereas the polarization drift is compressible. Hence, the 
inclusion of the polarization drift in the impurity model 
allows for genuinely compressible effects, such as those 
described in Sees. IIV( andlVl 



III. DIFFUSION OF IDEAL IMPURITIES 



B. Impurity passive-fluid model 

Impurities are here regarded as a passive species, that 
is, they are advected by the turbulence but do not have 
any influence on it. This is a reasonable approximation 
for the situation in which the density of impurity par- 
ticles is much lower than that of the bulk plasma. The 
modeling of the dynamics of impurities is consistent with 
that of the turbulence, in that it also follows from the 
drift-ordered fluid approach and the assumption of un- 
perturbed, uniform magnetic field. 

The impurities are characterized by a density and 
a flow velocity V6>. Just as in the case of the ions in 
the bulk plasma, the impurities are assumed cold, with 
their motion parallel to the magnetic field likewise being 
neglected. Therefore, the motion of the impurity fluid 
comprises only the electric drift v^; and the polarization 
drift Vp. Under the same normalization applied to the 



Theoretical studies of the diffusive effect of plasma 
edge turbulence on passively advected impurities have 
previously relied on analysis of the absolute dispersion 
of test particles The present approach is necessarily 
different, as the impurities are not modeled here as parti- 
cles but as a fluid. In this section we investigate the tur- 
bulent diffusion of ideal impurities by relative-diffusion 
analysis of the evolution of impurity puffs 

A. Relative difl'usion 

We consider the release of a puff of ideal impurities into 
the turbulent field, which is a saturated turbulent state 
described by the HW model. The amount of impurities 
released is 
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and the position c of the center of mass of the impurity 
puff at any time instant is given by 



evolution of Su can be characterized by effective diffusiv- 
ities Dj such that 



1 



xi9dx. 



We define the relative position vector = x — c, which 
indicates position in a reference frame moving with the 
center of mass of the puff. The average evolution of im- 
purity puffs can then be characterized by the following 
ensemble- averaged dispersion tensors: 



Sij 



1 



XiXjO dx 



Mi, 



{CiCj). 



i9dx 



Here, the components of vectors and tensors have been 
indicated by subindices and the angular brackets stand 
for the ensemble average, that is, the average over many 
releases. The absolute-dispersion tensor T^ij measures the 
dispersion with respect to the origin of the fixed frame, 
whereas the relative-dispersion tensor Sij measures the 
dispersion relative to the center of mass. The tensor Mij 
quantifies the "meandering" of the center of mass. The 
dispersion tensors are related by 

Yjij = Sij -\- Mij , 

which means that absolute dispersion consists of rela- 
tive dispersion and meandering of the center of mass. 
Relative-diffusion analysis focuses on the evolution of Sij , 
that is, on the spreading of the puff relative to its center 
of mass. Specifically, it is usually the diagonal elements 
Sii that are of interest once the reference axes are ade- 
quately selected, in our case those corresponding to the 
radial and poloidal directions. In view of the above rela- 
tion, the meandering of the center of mass is eliminated 
in relative-diffusion analysis, which thus leads to better 
statistics in comparison with those of absolute diffusion. 

The evolution of the dispersion tensors depends on the 
initial size and shape of the puffs, the statistical prop- 
erties of the turbulent velocity field, and the collisional 
diffusivity /j^q. Under the assumption that fig is small, 
the evolution of Sa for initially concentrated puffs fol- 
lows three stages: first, a phase dominated by collisional 
diffusion that lasts while the puff is much smaller than 
the turbulent structures and is characterized by 



dSu 
dt 



second, a transitional phase in which the size of the puff 
is comparable to that of the turbulent structures and Su 
typically evolve faster than linearly with time (see Ref.llSl 
for further discussions); and last, a phase dominated by 
turbulent diffusion that starts when the puff is large com- 
pared to the turbulent structures and during which the 



dSg 
dt 



2D,. 



(3) 



These diffusivities solely depend on the statistical prop- 
erties of the turbulent velocity field — more precisely, 
they are related to the mean squared velocities and the 
Lagrangian relative time scale. 

The last phase, dominated by turbulent diffusion, is 
most relevant for us, as it eventually dominates the trans- 
port of impurities. In fact, the effective diffusivities Di 
that characterize the phase of turbulent diffusion are 
equivalent in their definition to those employed in the ex- 
perimental modeling of anomalous impurity transport 
Because in this last stage the diagonal elements Su grow 
linearly according to Eq. , the effective diffusivities Di 
can easily be obtained from the asymptotic evolution of 
Su. This is the procedure used here for estimating the 
effective diffusivities of the turbulent field acting on ideal 
impurities. In principle, any initial puffs could be used for 
the determination of the diffusivities. From the discus- 
sion above, large initial puffs rapidly experience the final 
phase of turbulent diffusion, and thus provide the effec- 
tive diffusivities shortly after their release. However, care 
must be taken in numerical simulations that the puffs do 
not reach the boundaries of the simulation domain, as 
this would result in underestimates for the diffusivities. 
In order to determine the effective diffusivity in one di- 
rection, we take Gaussian stripes in the perpendicular 
direction as initial puffs and use periodic boundary con- 
ditions. In this way, the puffs are as large as possible in 
one direction while the conclusions of the simulations re- 
main valid. Moreover, because Gaussians are similarity 
solutions to the diffusion equation, the evolution of Su 
would be perfectly linear if the effect of turbulence were 
purely diffusive from the start. 



B. Simulation results 

The HW model for plasma edge turbulence, Eqs. (QJ, 
and the model for ideal impurities, Eq. Q with ( = 
0, were simultaneously solved numerically using the 
schemes described in the Appendix. The simulations 
were carried out on a doubly periodic square domain of 
side length 40 using 512 x 512 grid points. The viscosity 
/i^ and the collisional diffusivities /i^ and jj^q were set to 
the common value 0.02 in all the simulations presented 
in this paper. We considered four different values for the 
adiabaticity parameter C in the HW model: 0.5, 1, 2, 
and 4. As noted in the preceding section, low values of C 
correspond to hydrodynamic-like regimes, whereas high 
values correspond to nearly adiabatic regimes. In this re- 
spect, we regard C = 1 as a transitional value, representa- 
tive of an intermediate regime. An initial turbulent state 
was achieved by letting a uniform distribution of small- 
amplitude, random-phase modes evolve for a sufficiently 
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long time in order to reach a quasistationary saturated 
turbulent regime. Gaussian stripes having standard devi- 
ation 2 were used as initial impurity puffs. For each value 
of the adiabaticity parameter and each direction, radial 
and poloidal, 20 releases of such stripes were simulated, 
and the dispersion tensors were subsequently averaged. 

In Fig. n we present snapshots of the evolution of two 
perpendicular stripes of impurities released into the same 
saturated turbulent field with C = 1. The anisotropy in 
the diffusion of impurities is evident in the figure, as the 
poloidal spread of the puff in Fig. nje) is significantly 
larger than the corresponding radial spread in Fig.^^f). 
Hence, diffusion acts in this case faster in the poloidal 
direction than in the radial direction. The degree of 
anisotropy is in fact controlled by the adiabaticity pa- 
rameter, the turbulence ranging from isotropic in hydro- 
dynamic limit, C ^ 0, to strongly anisotropic in the adi- 
abatic limit, C oo. — 

The evolution of the diagonal elements Sa of the rela- 
tive dispersion tensor for the four values of the adiabatic- 
ity parameter is shown in Fig. [21 For the three smallest 
values of C, it is possible to distinguish the transitional 
and turbulent phases of growth in the plot for relative 
radial dispersion. In the plot corresponding to poloidal 
dispersion there is no clear distinction, as the evolution of 
5*22 very rapidly becomes approximately linear. For the 
side length of our simulation domain, dispersion can start 
becoming underestimated when Su reach values of about 
65. Nevertheless, for the three smallest values of C the 
linear trend in Su characterizing the phase of turbulent 
diffusion is established before such a point is reached. In 
contrast, the eventually linear regime is neither appar- 
ent in the radial nor the poloidal dispersion plots for the 
highest value of the adiabaticity parameter, C = 4. This 
different behavior is attributed to the weaker turbulence 
at higher values of C and the particular transport features 
of the adiabatic limit, the Hasegawa-Mima equation^ 

Estimates for the effective diffusivities Di of the turbu- 
lence for the three smallest values of C are obtained from 
linear fits to the linear portions of the curves in Fig. [21 
The estimated values for the radial, I^i, and poloidal, 1^2, 
effective diffusivities are presented in Tabled The values 
obtained here are qualitatively similar to those calculated 
in Ref. by means of a test-particle, absolute-diffusion 
approach. In particular, the decrease of Di with C as well 
as the rise of D2 with C for C > 1 are here reproduced. 
We note that the accuracy of the present approach could 
be improved by simulating the evolution of the impurity 
puffs for longer times on a larger domain. The number 
of puff releases could as well be increased in order to en- 
sure convergence to the ensemble averages. We also note 
that different values for ja^ and jUn were used in the lat- 
ter reference, with collisional diffusion of impurities being 
discarded because of their test-particle modeling. 




(e) (f) 



FIG. 1: (Color online) Evolution of Gaussian stripes of ideal 
impurities in a HW saturated turbulent state with C = 1: 
radial stripe at (a) t = 0, (c) t = 6, and (e) t = 12; and 
poloidal stripe at (b) t = 0, (d) t = 6, and (f) t = 12. 



TABLE I: Effective diffusivities in the radial direction Di 
and the poloidal direction D2 of HW saturated turbulence 
for different values of C. 



c 


Di 


D2 


0.5 


1.12 


2.32 


1 


0.65 


2.26 


2 


0.36 


2.63 



IV. CLUSTERING OF INERTIAL IMPURITIES 

The importance of inertia in the advection of particles 
by turbulent flows is well known in fluid dynamics. As 
in magnetized plasmas, the transport of inertial particles 
by fluids shows compressible features even when the ad- 
vecting flow is incompressible. Particle inertia has been 
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FIG. 2: Evolution of the relative radial dispersion Sn and the 
relative poloidal dispersion S22 of ideal impurity Gaussian 
stripes released in HW saturated turbulence with C = 0.5 
(solid), C = 1 (dashed), C = 2 (dotted), and C = 4 (dash- 
dotted). 

shown to result in clustering in vortical structures, and its 
role in phenomena such as planetary formation has been 
suggested. While the dynamics of impurity particles in 
magnetized plasmas are not governed by the same forces 
as in fluids, inertia enters the dynamics of impurities in 
magnetized plasmas in a similar way as the Coriolis force 
enters those of heavy impurities in rotating fluids. Hence, 
corresponding clustering effects for inertial impurities in 
magnetized plasmas should be expected. In this section 
we discover and explain the correlation between impurity 
density and vorticity that results from the compressibility 
introduced by the polarization drift in the case of inertial 
impurities. In order to most clearly reveal this effect we 
consider a homogeneous initial distribution of impurities, 
since such a distribution would remain unchanged in the 
incompressible, that is, ideal-impurity case. 



A. Simulation results 

Simulations of the HW model for plasma edge turbu- 
lence, Eqs. and the impurity model, Eq. (|2j, were 
again performed on a doubly periodic square domain of 
side length 40 using 512 x 512 grid points. The adia- 
baticity parameter C in the HW model was set to the 



transitional value 1, while we contemplated various val- 
ues for the parameter ranging from —0.01 to 0.05. The 
impurity density field was initialized with the constant 
value ^0 = 1, and the initial turbulent fields were taken 
from a HW saturated turbulent state. 

In Fig.iniwe present side by side snapshots of the evo- 
lution of the vorticity uj of the electric drift velocity and 
the impurity density field for the case ( = 0.01. Merely 
5 time units after the simulation is started, there is al- 
ready indication of a correlation between vorticity and 
impurity density. At t = 50, there is visually little differ- 
ence between the two fields. Figure0]shows a scatter plot 
of vorticity and relative impurity density 6/60 dit t = 100 
for three different values of The plot clearly suggests 
an approximate linear relation given by 

e/Oo ^ 1 + i^cj, 

with the regression coefficient K depending on The 
least-squares estimates of the coefficient K are plotted 
as a function of the parameter ( in Fig. [5J This figure in 
turn indicates a linear relation between K and ( in the 
form K ^ 0.82^. In summary, we find that the density 
of inertial impurities eventually becomes linearly corre- 
lated with the vorticity, the regression coefficient being 
proportional to the parameter 

B. Turbulent mixing and clustering 

The behavior just described can easily be explained 
using an argument of turbulent mixing, or turbulent 
equipartition^2i2ii2S Performing basic algebraic manip- 
ulations on Eq. ^ while accounting for the particular 
definition of the Lagrangian derivative in Eq. (|ld|) . 
we may rewrite the impurity model equation in the form 

D,(ln^ - Coo) = CVi_ \nO • D,V^(/9 + ^ Vi^, (4) 

u 

where the left-hand side describes advection by the elec- 
tric drift and compression by the polarization drift, while 
the right hand-side describes advection by the polar- 
ization drift as well as collisional diffusion. We will 
now argue that the right-hand side of this equation is 
small, which implies that \nO — (co is approximately a 
Lagrangian invariant, or in other words, that \nO — (co 
is nearly conserved along the trajectories defined by the 
electric drift velocity. We split the impurity density 
into its mean value Oq and the fluctuations 61 arising from 
compressibility. Based on the computational results, we 
further postulate that Oi/Oq = 0(C)- Selectively apply- 
ing the split = Oq -\- Oi to Eq. (Q, we obtain 

Dtilne- Cc^) = CVi ln(^l + • D^V^^ + ^ Vl^i. 

Assuming that the normalizations used in the impurity 
model adequately represent the scales of fluctuations, the 





















\ 





FIG. 3: (Color online) Evolution of the vorticity u and the 
density of inert ial impurities (C = 0.01) in a HW saturated 
turbulent state with C = 1: (a) a; and (b) ^ at t = 0, (c) oo 
and (d) ^ at t = 5, and (e) oo and (f) ^ at t = 50. Green corre- 
sponds to mean values, while red and blue represent positive 
and negative fluctuations respectively. 



right-hand side of this last equation is of order + (fig^ 
whereas the fluctuating component of the quantity inside 
the Lagrangian derivative on left-hand side is of order 
Hence, \nO — (co varies slowly along the electric drift tra- 
jectories, thus constituting an approximate Lagrangian 
invariant. 

Turbulent mixing homogenizes Lagrangian invariants, 
and thus tends to drive systems to states of so-called tur- 
bulent equipartition. The application of this argument to 
the approximate invariant here yields 

\nO — (co ^ const. 

Because the amount of impurities is conserved and it is 
assumed that 6i/0o = 0{C), the previous expression can 
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FIG. 4: (Color online) Scatter plot of vorticity co and relative 
density O/Oq of inertial impurities in a HW saturated turbu- 
lent state with C = 1 for = 0.05 (red/steepest), C = 0.01 
(green), and C = 0.002 (blue/flattest). 
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FIG. 5: Least-squares estimates of the coefficient X as a func- 
tion of C for the impurity-density-vorticity linear regression 
O/Oo ^ l-\-K(jj in a HW saturated turbulent state with C = 1. 
The fitted line (dashed) corresponds to K = 0.82^. 



be recast into the simpler form 

O/Oo ^ 1 + c^. 

The present analysis thus predicts a linear relation be- 
tween the impurity density and the vorticity, the regres- 
sion coefficient being the parameter (. Hence, the the- 
oretical argument explains the effect discovered in the 
simulations up to a slight mismatch in the regression co- 
efficient. We emphasize that the particular features of 
the HW turbulence model were not used in the previous 
argument, the conclusions relying solely on the modeling 
of the dynamics of impurities and the argument of tur- 
bulent mixing. It follows that the clustering of inertial 
impurities in vortical structures is a generic effect, in- 
dependent of the underlying instability mechanisms and 
the specific characteristics of the turbulence. 

The fact that the density of inertial impurities becomes 
directly correlated with the vorticity implies that impu- 
rities of positive charge cluster in positive vortices, the 
opposite taking place for negatively charged impurities. 
As previously introduced, this effect is similar to the ag- 
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gregation of heavy particles in anticyclonic vortices that 
can take place in rotating fluids as a result of the Coriolis 
forced In contrast with compression by the polarization 
drift, which yields D^ln^ ~ C^t^^ compression by the 
Coriolis force enters the dynamics of heavy impurities in 
the form D^ln^ ~ —2Qij^ with Q being the overall rota- 
tion. Hence, aggregation in the rotating-fluid case does 
not necessarily lead to a linear relation between impu- 
rity density and vorticity, but can go on further to form 
highly concentrated clusters in the cores of anticyclonic 
vortices. For this reason, it has for example been sug- 
gested as a mechanism, in combination with gravitation, 
for planetary formation in rotating astrophysical disks. 



V. PINCH OF INERTIAL IMPURITIES 

In this section we finally investigate the possible role 
of inertia in the radially inward pinch that concentrates 
impurities in the center of the plasma column in magnetic 
confinement devices. As a measure of the overall drift 
of impurities, we consider the global impurity fiux 



J Owe dx + j 



Ovp dx. 



Here, we have distinguished the components respectively 
resulting from advection by the electric drift and by the 
polarization drift. In terms of these global fiuxes, we 
define the impurity drift velocity \e and its components 
Vf and by means of 



-ItE 



-i-rP 



where once again Q = J dx. For a homogeneous distri- 
bution of impurities in a periodic domain, both Tq and 
Tg are zero. Hence, ideal impurities cannot experience a 
net drift when their initial distribution is homogeneous, 
since such a distribution remains unchanged in the ideal 
case. In contrast, inert ial impurities may experience a 
drift even if their initial distribution is uniform, given 
that the latter is altered by compressible effects, as shown 
in the preceding section. 



A. Simulation results 

We begin by monitoring the radial drift velocity \e ' x 
in the simulations referred to in Sec. II VI These were 
initialized with homogeneous distributions of inertial im- 
purities in a HW saturated turbulent state with C = 1. 
In Fig. we present the evolution of the radial drift ve- 
locity for four different values of (. In every case, the 
evolution of the radial drift comprises a strong transient 
burst, after which the drift enters a saturated quasista- 
tionary regime. Moreover, the radial drift velocity is seen 
to have a definite sign opposite to that of that is to say, 
opposite to the type of charge of the impurity particles. It 
follows that positively charged impurities are subject to 
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FIG. 6: Evolution of the radial drift velocity • x of inertial 
impurities in a HW saturated turbulent state with C = 1 for 
( = -0.01 (solid), C = 0.005 (dashed), C = 0.02 (dotted), and 
C = 0.05 (dash-dotted). 
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FIG. 7: Time- averaged radial drift velocity • x of inertial 
impurities as a function of C in a HW saturated turbulent 
with C = 1. The fitted line (dashed) corresponds to V6» • x = 
-0.090C. 



a continuous radially inward drift, as observed in exper- 
imental investigations of impurity transport. ^'^ In Fig.Q 
we show the time-averaged radial drift velocities \e ' x 
in the saturated regimes, computed using the values of 
the drifts between t = 25 and t = 150. The variation of 
the time-averaged radial drift with ( is remarkably linear, 
well fitted by • x = — 0.090C. Therefore, the present 
radial drift scales linearly with the mass-charge ratio of 
the impurity particles. 

The reason for this peculiar pinch effect is not clear 
at the outset. An analysis of the components and 
\q of the drift velocity in the simulations reveals that 
the contribution of the electric drift is by far dominant. 
Hence, while the net drift arises as a result of the polar- 
ization drift, the main role of the latter is not to globally 
advect impurities, but rather to distribute them in such 
a way that they experience net transport by the electric 
drift. From a naive point of view, it could be thought 
that the correlation between impurity density and vor- 
ticity discovered in Sec. IIVI is responsible for the radial 
drift. Indeed, transport of trapped particles by coher- 



9 



ent vortices is for instance known to play an important 
role in transport processes in rotating fluids. However, 
because an incompressible flow in a 2D periodic domain 
cannot globally transport its vorticity, a purely linear re- 
lation between impurity density and vorticity would yield 
no drift at all. Thus, the explanation for the radial drift 
of inertial impurities requires deeper thought. 



B. Net transport and pinch 

Intrigued by the initial burst and the subsequent sat- 
uration of the radial impurity drift velocity, we analyze 
the evolution of the global impurity flux. In particular, 
we focus on the Tq component, given that the net trans- 
port of impurities is dominated by the electric drift. In a 
periodic domain, the evolution of Tq under the impurity 
model in Eq. (|2|) is governed by 



dt 



I 



Ai = / ODt^rE dx. 



where 



C6>V£;DtCjdx, 

A3 = J C{'^±0'DtV^^)vE dx, 

A4 = y /i0V£;Vil9dx. 

These four contributions correspond respectively to evo- 
lution of the electric drift velocity, compression by the po- 
larization drift, advection by the polarization drift, and 
collisional diffusion of impurities. In order to reveal the 
magnitude of the different terms, we again split the im- 
purity density into its mean value ^0 and the fluctuations 
^1, keeping in mind that Oi/Oq = 0{Q. Selectively ap- 
plying the split = Oq -\- Oi to the above definitions we 
obtain 



Ai 
A2 
A3 
A4 



J ^iD,V£;dx=(9(C), 

I aV^Oi . D, V^(^)v^ dx = O(C'), 
J ^o^eVIOi dx = (9(C/i^). 



According to these estimates, Ai and A2 dominate the 
dynamics of the global impurity flux. In the case of uni- 
form impurity density, A2 is the only nonvanishing con- 
tribution to dr^/dt, being therefore responsible for the 
initial development of drift velocities. 



We will now resort to the particular HW turbulence 
model to show that A2 constitutes a source for global 
radial flux in the present case. Substitution of the HW 
vorticity equation in the above expression for A2 yields 

A2 = j CO[C{^-n)^ ii^V\uj\wEd-^ 

= -CC^o j nwE dx + (9(C' + CMo.). 

Projecting this equation onto the radial direction we find 

A2 • X = -CC6>o j ndyLpd^ + 0{C^ + Cm^) 

Here, we have identified the sole source of energy and 
generalized enstrophy within the HW model. Since is 
in practice positive definite (see, e.g.. Table I in Ref . Illl) . 
it follows that A2 • x has a definite sign opposite to that 
of C, thus acting as a source for the global radial flux 
of inertial impurities. As explained before, A2 is initially 
the only nonzero contribution to dPf /dt, thus explaining 
the transient bursts observed in Fig. El More generally, 
this term sustains the radial transport of impurities in 
the turbulent states described by the HW model. 

From the above estimates, we are led to expect that 
Ai is responsible for the saturation of the radial impu- 
rity drift seen in Fig. |H1 The reason for this saturation 
becomes clear when we account for the correlation that 
develops between impurity density and vorticity. The 
argument of turbulent equipartition used in Sec. II VI pre- 
dicted an eventual relation Oi/Oq ^ between the impu- 
rity density fluctuations and the vorticity. Considering 
this relation and the fact that vorticity is not globally 
transported by the electric drift we obtain 



Ai 



/ 



COoujDtVE dx 



J C^0V£;D 



tuj dx ! 



whenever turbulent equipartition holds true. Conse- 
quently, when such a state is reached, the two terms 
dominating the evolution of the radial drift approxi- 
mately cancel, and the radial drift enters a quasista- 
tionary regime. The time it takes for this regime to be 
reached is independent of as the development of corre- 
lation is exclusively dependent on the turbulent mixing 
by the electric drift. In fact, the evolution of Oi is approx- 
imately linear with ( within short time intervals, since in 
line with the discussion in Sec. |^ it follows that 

Given that solely the density fluctuations yield net trans- 
port of impurities, the evolution of the impurity drift is 
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likewise approximately linear with This explains the 
apparent linear scaling of the curves in Fig. El and the 
linearity of the time-averaged radial drift in Fig. [3 

While the arguments presented here do not yield a pre- 
diction for the quasistationary level of the radial impurity 
drift, they do identify the mechanism for its onset and 
likewise explain its scaling with the parameter (. We 
emphasize that the former mechanism relies on the spe- 
cific properties of the turbulence described by the HW 
model. Nevertheless, an analogue pinch effect would ap- 
pear for other driving instabilities provided that the term 
A2 remains a source for a definite impurity drift. 



charge ratio of the impurity particles and is inward for 
positively charged impurities. This anomalous pinch re- 
lies on the particular features of the turbulence described 
by the HW model. In contrast with other known turbu- 
lent pinches,^ the present pinch is neither caused by a 
specific magnetic field geometry nor by temperature gra- 
dients, given that our models assume a uniform magnetic 
field and cold impurities. Consequently, our finding of a 
net radial drift due to impurity-particle inertia consti- 
tutes a new contribution towards the understanding of 
the anomalous pinch of impurities in magnetic confine- 
ment devices. 



VI. CONCLUSION 

We have studied the transport of impurity particles in 
plasma edge turbulence, basing our investigations on the 
paradigmatic Hasegawa-Wakatani (HW) model for resis- 
tive drift-wave turbulence and a consistent passive-fluid 
model for the impurities. The latter model takes into 
account impurity-particle inertia, which enters the im- 
purity flow velocity in the form of the polarization drift. 
Inertia has indeed been found to play a significant role in 
the turbulent transport of impurities, since it gives rise to 
subtle yet qualitatively important compressible effects. 

The model equations have been investigated both the- 
oretically and computationally. The numerical simula- 
tions were performed using the semi-Lagrangian (SL) 
pseudospectral code briefly described in the Appendix. 
Although the emphasis of this paper is not on the nu- 
merical methods employed, we point out the suitability 
of SL schemes for the advection-dominated problems that 
appear in plasma turbulence. In particular, we highlight 
the superior stability of the SL semi-implicit scheme used 
here for the HW model, which is known to constitute a 
numerically challenging nonlinear problem. 

In the case of ideal impurities, we have focused on their 
turbulent diffusion by the electric drift. This effect has 
been demonstrated by relative-diffusion analysis of the 
evolution of impurity puffs. The effective diffusivities 
of the turbulent advecting field have been calculated for 
several values of the adiabaticity parameter in the HW 
model. The resulting turbulent diffusivities are in quali- 
tative agreement with those obtained in an earlier study 
based on an absolute-diffusion, test-particle approach. 

In the case of inert ial impurities, we have discovered 
that their density eventually correlates with the vorticity 
of the electric drift velocity. This clustering effect scales 
linearly with the mass-charge ratio of the impurity par- 
ticles and results from compression by the polarization 
drift in combination with turbulent mixing. The cluster- 
ing of impurities in vortices of definite sign is a generic 
effect in magnet ized-plasma turbulence, in the sense that 
it is independent of the specific characteristics of the tur- 
bulence. 

Finally, we have also found an overall radial drift of 
inertial impurities that scales linearly with the mass- 
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APPENDIX: NUMERICAL SCHEMES 

The numerical simulations presented in this paper 
were performed using a semi-Lagrangian (SL) pseudo- 
spectral code developed by the first author^ In essence, 
SL schemes combine the very stable Lagrangian treat- 
ment of advection with the convenience of a fixed spa- 
tial discretization^ii2& This is achieved by integrating the 
partial differential equations in question along a different 
family of advective trajectories in each time step, pre- 
cisely those that by the end of the time step traverse the 
fixed points on which the spatial discretization is based. 
We now outline how the aforementioned code combines 
SL method and the pseudospectral discretization to nu- 
merically solve the HW model and the impurity passive- 
fluid model described in Sec. ^ 

We discretize time using a constant time step At, 
which yields the sequence {tm} of time instants related 
by tm-\-i = tm-\-^t. Likewise, we use the Fourier pseudo- 
spectral method to discretize the doubly periodic rect- 
angular simulation domain. This results in a regular 
grid {x^} of collocation points. In our application of 
the SL method we consider the (incompressible) electric 
drift velocity v^; = z X V^(/? as the sole advecting field. 
Hence, the advective trajectory X(x^, tm+i; t) that passes 
through the grid point at time t^+i is given by 

^X(Xi,t^+i;t) = V£;(X(Xi,tm+i;t),t), 
X(Xi, trn+i; ^m+l) = ^i- 

At each time step, we firstly determine the so-called up- 
stream points {X(x^, t^+i; t^) = x^ — a^}. These are 
numerically calculated using the following second-order 
scheme, which combines the implicit midpoint rule and 
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linear extrapolation of the advecting fieldjSSiSi 



AtvJ;(> 



^m+1/2) 



vJ;(x,t^+l/2) = fv£;(x,t^) - |v£;(x,t^_i), 

where = + 1 At. The implicit equations for the 

displacements {a^} are solved using Newton's method, 
which converges in very few iterations. 

Subsequently, the model equations are integrated from 
tm to along the advective trajectories. For the HW 
model, Eq. (^P), we use the following unconditionally sta- 
ble second-order semi-implicit scheme: 



- n)™+i + - n) 



,m+l 



Here the superscripts indicate time, and tildes denote 
evaluation at the upstream points. Hence, for a func- 
tion g(x, t) we have g^(xi) = g(xi — a^, tm). The above 
equations supplemented with uo^^^ = \/']_(p^~^^ are eas- 
ily solved within the Fourier discretization. 



For the impurity passive-fluid model, Eq. ((^J, we em- 
ploy the following second-order implicit-explicit scheme: 



gm+l nm 



In our code we further approximate the polarization drift 
Vp = — CDtVx^ by means of 



C(VI^^ - V±ip'^^^)/At + (9(At). 



While in theory this approximation reduces the accuracy 
of the scheme to flrst order, its effect in practice is very 
limited given the smallness of the polarization drift. 

We lastly note that the combination of spectral dis- 
cretizations with SL schemes is problematic^ The diffi- 
culty lies in the required evaluations at midpoints and up- 
stream points, which are generally not equispaced. Per- 
forming such evaluations exactly, using the spectral rep- 
resentation, is in most cases computationally infeasible. 
In our code this problem is circumvented using high-order 
periodic spline interpolation on a refined grid»22iS2. 
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